function [yobsCounter,statesCounter]=kfilterNSplitsCounterfactual(parvecNew,funcmod,Y,solveopt,...
    addsol,etaOrig,aZero)
% ====================================================================
%% kfilterNSplitsCounterfactual
% 
% Given the original set of shocks inferred with the original parameter, 
% solve the model with a new parameter parvecNew, feed in the original set
% of shocks, and obtain the counterfactual evolution of the states and the
% observables, with the constant added back in. 
% 
% Must pass is *aZero* the original smooth guess of the states, such that 
% if parvecNew the original parameter vector, this will recover the
% observables. 
% ====================================================================

%% 1. Model solution
% Matrices of second sample stored in structure second 
[G,R,C,eu,SDX,Z]=feval(funcmod,parvecNew,solveopt,addsol);
if isequal(eu,[1;1])==0 
    error('Model is not determinate') 
end 

%% 2. Forward filter 
% 3.1. Dimensions and storage 
[T nobs]=size(Y); 
ns = size(G,1); 
nx = size(SDX,1); 

%% 3. Smoother with original shocks 
% Check that Smooth States are identical if using disturbance smoother (above) vs. state smoother (below) 
% and if can also recover the observables
% State at time zero give by GG1*azero+Pzero*rstar;
[yobsCounter,statesCounter]=kfilterRegSplitSimulation(aZero,etaOrig'); 


%% Subroutine kfilterRegSplitSimulation allows to simulate the model  
% Inputs 
% sInitial: [ns 1] Initial State 
% innovMat: [nx T] matrix of innovations 
% By being a sub-routine it has access to all variables defined above 
% Be-careful not to use an index (ii,jj) used above or to repeat variable
% names
    function [ySim,sSim]=kfilterRegSplitSimulation(sInitial,innovMat) 
        if ~isequal(size(innovMat),[nx T]) 
            error('Input innovMat must be [nx T]') 
        end 
        sSim=zeros(ns,T);
        ySim=zeros([nobs T]);
        sSim(:,1)=G(:,:,addsol.tauVec(1))*sInitial + ...
            R(:,:,addsol.tauVec(1))*innovMat(:,1);
        ySim(:,1)=Z(:,:,addsol.tauVec(1))*(sSim(:,1)+C(:,addsol.tauVec(1)));
        for kk=2:T;
            sSim(:,kk)=G(:,:,addsol.tauVec(kk))*sSim(:,kk-1)+...
                R(:,:,addsol.tauVec(kk))*innovMat(:,kk);
            ySim(:,kk)=Z(:,:,addsol.tauVec(kk))*(sSim(:,kk)+C(:,addsol.tauVec(kk)));
        end
        ySim=ySim'; 
        sSim=sSim'; 
    end
%% End of File
end 